AD-A277  648 


NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CAUFORMA 


THESIS 


DESIGN  AND  MONTE  CARLO  ANALYSIS  OF  AN 
ON»MnED  AERIAL  VEHICLE 


Josei^i  PbuI  Foidham 
December  1993 


ris  Advisor. 


Isaac  Kaminer 


^^roved  for  wblic  release;  dtstribution  is  unlimited 


cMi  94-09959 

■nniin 


DTTC  QTTP^Vmr  W^WnTED  1 


94  4  1  O40 


REPORT  DOCUMENTATION  PAGE 

form  Afiprovad 

0MB  No.  0704-dtBB 

PuWfC  ft9ontt*9  bur^  for  tbo  coUoctioj^  information  it  ottimattO  to  avora^o  t  nowr  otr  raiponta.  mdudino  tfit  tima  for  ravitwiwo  inttrwcttont.  itarchino  eaittmo  Oata  toorctt 
^thtrmo  and  maintaining  tha  Oata  nacOad.  and  comoiating  and  ra^iawing  tha  coiiactton  of  information  Sand  commanti  ragarding  tmt  iMrtton  attimata  or  any  othor  atc^  of  ttwt 
CQllactiQn  of  information,  indudmg togMtiont  for  raducing  enn  bordan.  lo  Watnmgton  Haadouartan  Sarvicat.  Oiraaorata  for  information  Ooaratiom  and  Raot^  MM  Jaffarton 
Omni  Highway.  Suita  1204.  Arimgten.  VA  22202*4302.  and  to  tna  Offica  of  Managamant  and  tudgat.  oapanwork  Radoction  Proiact(07044)it0).  dfaihington.  DC  20S0I. 

1.  AeiNCV  USt  ONLY  a««vt  bisnk)  2.  RfAORT  DATE  3.  RCAORT  TYRE  AND  OATES  COVERED 

December  1993  Master's  Thesis  4/93-10/93 

4.  TITLE  AND  SUITITLE 

Design  and  Monte  Carlo  Analysis  of  an  Unmanned  Aerial 
Vehicle 

S.  FUNDING  NUNWERS 

S.  AUTHOR(S) 

Joseph  Paul  Fordham 

7.  RERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS(ES) 

Naval  Postgraduate  School 

Monterey,  California  93943-5000 

B.  RERFORMING  ORGANIZATION 
RERORT  NUMBER 

9.  SRONSORING/ MONITORING  AGENCY  NAME(S)  AND  AOORESS(ES) 

10.  SRONSORING /MONITORING 

AGENCY  RERORT  NUMBER 

11.  SURRUMENTARY  NOTES 

12s.  DISTRIBUTION /AVAILARILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words) 

In  the  last  several  years,  software  innovations  and 
the  increasing  speed  and  availedsility  of  microcomputers 
and  workstations  have  made  the  dynamic  simulation  of  complex 
systems  more  practical.  One  such  system,  a  short-range  Unmanned 

Aerial  Vehicle  called  Bluebird,  was  previously  modeled  on  Simulink, 
a  commercial  software  package.  The  high  fidelity  model  includes  six 
degree  of  freedom  nonlinear  equations  of  motion  with  onboard  sensors 
and  a  Global  Positioning  System  and  inertial  navigation  system. 

Because  of  interest  expressed  by  the  Unmanned  Aerial  Vehicle 

Joint  Program  Office  in  how  accurately  a  UAV  could  identify  a  target's 
geographical  coordinates,  the  Bluebird  model,  with  an  added  guidance 
and  control  system,  was  evaluated  as  to  its  navigational  and 
attitudinal  accuracy  in  a  dynamic  simulation  using  Monte  Carlo 
techniques.  Because  of  the  modular  nature  of  the  simulation,  future 
evaluations  of  manned  or  unmanned  aircraft  and  avionics  will  involve 
relatively  uncomplicated  changes  to  the  existing  model. 

14.  SUBIEO  TERMS 

IS.  NUMBER  OF  RAGES 

89 

IB.  RRKE  CODE 

17.  SECURITY  CUSSIRKATION 
OR  RERORT 

UNCLASSIFIED 

IB.  SECURITY  CLASSIFICATION 
OF  THIS  RAGE 

UNCLASSIFIED 

19.  SECURITY  CUSSIRKATION 
OF  ABSTRAO 

UNCLASSIFIED 

20.  LIMITATION  OF  ABSTRACT 

UL 

NSN  754(M)1-280-S$00  I  Standard  Form  298  (R«v  2-89) 


PrCKfibcd  by  ANil  Sid  Ot-^t 
29a-l02 


Apfvoved  for  public  release;  discribution  is  unlimited. 

DESIGN  AND  MONTE  CARLO  ANALYSIS  OF  AN  UNMANNED  AERIAL  VEHICLE 


Author 


Approved  by; 


by 

Joseph  P.  Fordham 
Lieutenant.  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1985 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  AERONAUTICAL  ENGINEERING 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 


Michael  K.  Shields,  Second  Reader 


Daniel  J.  CollinM;  Chaimian 
Department  of  Aeronautics  and  Astronautics 


ii 


ABSTRACT 


In  the  last  several  years,  software  innovations  and  the  increasing  speed 
and  availability  of  microcomputers  and  woiicstadons  have  made  the  dynamic  simu¬ 
lation  of  complex  systems  more  practical.  One  such  system,  a  short-range 
Unmanned  Aerial  Vehicle  called  Bluebird,  was  previously  modeled  on  Simulink,  a 
commercial  software  package.  The  high  fidelity  model  includes  six  degree  of  free¬ 
dom  nonlinear  equations  of  motion  with  onboard  sensors  and  a  Global  Positioning 
System  and  inertial  navigation  system. 

Because  of  interest  expressed  by  the  Unmanned  Aerial  Vehicle  Joint  Pro¬ 
gram  Office  in  how  accurately  a  UAV  could  identify  a  target's  geogr^hical  coordi¬ 
nates,  the  Bluebird  model,  with  an  added  guidance  and  control  system,  was  eval¬ 
uated  as  to  its  navigational  and  atdtudinal  accuracy  in  a  dynamic  simulation  using 
Monte  Carlo  techniques.  Because  of  the  modular  nature  of  the  simulation,  future 
evaluations  of  manned  or  unmanned  aircraft  and  avionics  will  involve  relatively 
uncomplicated  changes  to  the  existing  model. 
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I.  INTRODUCTION 


In  the  last  several  years,  the  increasing  speed  and  availability  of  microcomput¬ 
ers  and  workstations  has  allowed  complex  simulation  of  dynamic  systems  which 
previously  could  only  be  practically  evaluated  at  static  points  in  time.  In  order  to 
rigorously  evaluate  the  performance  of  a  nonlinear  multiple  input,  multiple  output 
(MIMO)  modeled  air  vehicle,  such  as  an  unmanned  aerial  vehicle  (UAV),  it  is  nec¬ 
essary  to  dynamically  simulate  it.  This  is  due  to  the  complex  way  in  which  errors 
are  fnopagated  throughout  a  nonlinear  closed  loop  system.  The  most  common  way 
to  dynamically  evaluate  errors  of  a  system  subject  to  noise  is  using  the  Monte  Carlo 
simulation  method,  which  is  essentially  the  repetition  of  the  same  experiment  over 
many  times. 

To  sui^rt  the  ongoing  Naval  Postgraduate  School  UAV  effort,  a  high  fidelity 
UAV  computer  model,  incorporating  aircraft  equations  of  motion  and  sensors, 
guidance,  navigation,  and  control  models,  was  constructed  utilizing  the  thesis  work 
of  two  former  Naval  Postgraduate  School  students.  The  primary  goal  was  an  anal¬ 
ysis  of  how  well  a  UAV  could  designate  a  fixed  ground  target  considering  the  reli¬ 
ability  of  its  onboard  sensors.  To  achieve  this,  the  model,  with  a  simulated  FLIR 
(Forward  Looking  Infrared  Receiver)  or  low  light  television  payload  attached,  was 
run  over  a  fixed  track  on  a  reconnaissance  mission.  For  the  Monte  Carlo  simula¬ 
tion  and  analysis,  sensor  error  data  was  taken  from  a  Naval  Air  Development  Cen¬ 
ter  (NADC)  report,  "Unmanned  Aerial  Vehicle  Flight  Management  System  Error 
Artalysis"  [Ref.  1]  and  added  to  the  UAV  model.  The  Monte  Carlo  analysis  per- 
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fonned  on  this  model  was  intended  to  complement  [Ref.  1],  which  analyzed  and 
compared  targeting  errors  for  several  distinct  UAV  architectures,  including  close- 
range,  short-range,  and  medium-range  vehicles  performing  a  similar  mission.  This 
thesis  is  an  extension  of  the  UAV  Error  Analysis  Report  in  that  it  considers  how 
sensor  and  measurement  errors  {»x>pagate  through  a  UAV  guidance,  navigation,  and 
control  system  as  it  directs  the  vehicle  along  a  given  trajectory.  This  report  shows 
the  importance  of  considering  the  feedback  nature  of  a  controlled  system  in  the 
error  analysis  of  its  performance. 

Consider  the  system  shown  in  Figure  1.  Here  the  aircraft  block  contains  actua¬ 
tors  and  sensors  which  are  subject  to  measurement  and  position  errors.  Errors  from 
one  sensor  used  by  the  control  system  for  feedback  propagate  into  all  feedback 
channels  in  a  complex  nonlinear  manner.  For  example,  a  change  in  roll  rate  may 
cause  a  false  measurement  in  the  pitch  rate  gyro.  If  this  gyro  is  being  used  to  help 
determine  aircraft  attitude,  this  false  measurement  will  cause  an  inaccuracy  in  mea¬ 
sured  attitude.  To  restore  proper  attitude,  an  actual  pitch  rate  will  be  created,  which 
may  cause  a  false  measurement  of  roll  and  yaw  rates  and  correspondingly  of  roll 
and  yaw  angles.  It  is  due  to  the  nonlinear,  feedback  nature  of  these  phenomena  dial 
simple  statistical  methods  applied  to  a  static  model  of  a  vehicle  cannot  be  used  to 
accurately  evaluate  how  well  a  sensor  suite  measures  aircraft  flight  parameters.  In 
the  analysis  presented  here,  consideration  was  also  given  to  sensor  dynamics  (how 
sensitive  sensors  are  to  a  rapidly  changing  input).  Angular  sensors  on  the  aircraft 
and  lookdown  and  azimuth  errors  for  the  FLIR/  Low  Light  Television  sensor  were 
those  specified  for  the  MI  AG  architecture  in  [Ref.  1].  To  complete  the  analysis, 
Monte  Carlo  simulation  was  implemented  with  the  NFS  UAV  model  by  repeatedly 
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running  the  model  over  a  fixed  track  and  collecting  and  analyzing  the  onboard  sen¬ 
sor  data.  The  targeting  errors  obtained  given  the  advertised  sensor  accuracies  were 
then  compared  with  the  taigeting  errors  determined  with  Monte  Carlo  simulation. 

To  accomplish  the  error  analysis,  the  UAV  model  developed  by  two  prior  thesis 
students  had  to  be  completed.  Figure  1  shows  how  their  work  fits  into  the  complete 
model.  The  model  consisted  of  aircraft  nonlinear  equations  of  motion  and  sensor 
models  developed  by  Kuechenmeister  [Ref.  2]  and  a  high  fidelity  Differential  Glo¬ 
bal  Positioning  System  and  inertial  navigation  system  model  created  by  Marquis 
[Ref.  3].  A  linear  quadratic  regulator,  whose  design  and  implementation  is  dis¬ 
cussed  in  Chapter  HI,  was  constructed  to  stabilize  the  model  and  provide  a  means 
to  control  it  in  heading,  airspeed,  and  attitude  in  steady  state.  To  steer  the  model 
along  a  fixed  course,  a  waypoint  guidance  feature  was  also  developed.  While  the 
entire  model  is  not  representative  of  a  particular  UAV,  the  way  that  its  guidance, 
navigation,  and  control  are  interrelated  in  their  operation  is  indicative  of  how  a 
UAV  equipped  with  an  integrated  sensor  package,  like  the  new  MIAG  (modular 
integrated  avionics  group)  concept,  can  be  sensitive  to  sensor  errors. 

The  next  chapter  outlines  the  components  of  the  computer  model.  Chapter  III 
covers  in  greater  depth  the  coontributions  to  the  model  made  by  the  author,  includ¬ 
ing  controller  and  waypoint  guidance  design  and  payload  modeling.  Chapter  FV 
explains  the  Monte  Carlo  analysis  used  in  analyzing  the  model’s  taigeting  errors. 
Chapter  V  discusses  how  the  simulation  was  performed  and  Chapters  VI  and  VII 
include  results  of  the  analysis  and  conclusions. 
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Figure  1.  System  Overview 


II.  THE  UAV  MODEL 

The  complete  model  used  in  this  thesis  is  shown  in  Figure  2.  It  is  composed  of 
several  distinct  subsystems,  some  of  which  were  constructed  as  parts  of  previous 
theses  by  Naval  Postgraduate  School  Students.  Kuechenmeister  was  responsible  in 
a  large  part  for  the  aircraft  equations  of  motion  and  the  sensor  models  [kcr.  2]. 
Marquis  constructed  the  six  satellite  differential  Global  Positioning  System  (GPS) 
model  used  in  the  simulation  as  well  as  a  Kalman  Filter  integrating  the  GPS  model 
with  an  inertial  navigation  system  [Ref.  3].  The  model,  which  was  constructed 
using  Simulink,  a  program  commercially  available  from  Mathworks,  Inc.,  is 
designed  to  run  on  either  a  Sun  workstation  or  a  personal  computer  miming  Matlab 
for  Windows  with  SimuUnk.  Because  of  the  complexity  of  the  simulation,  it  should 
be  mn  on  the  fastest  platform  available.  As  a  benchmark,  the  entire  simulation  mn- 
ning  on  a  Sun  SparclO  workstation  with  64  megabytes  of  RAM  rans  at  a  rate  of 
about  one  second  of  real  time  for  every  five  minutes  of  simulation  time. 

A.  EQUATIONS  OF  MOTION  AND  SENSOR  MODELS 

The  UAV  used  in  the  simulation  is  called  Bluebird.  The  actual  aircraft  is  a  high- 
wing  monoplane  with  a  wingspan  of  12.4  feet  and  a  weight  of  55  pounds.  The  air¬ 
craft  is  controlled  using  servos  which  actuate  conventional  elevator,  aileron,  and 
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Figure!.  System  Diagram 

rodder  surfaces  and  a  throttle  controlling  a  half-horsepower  piston  engine.  Motive 
force  is  provided  by  a  single  nose-mounted  tractor  propeller.  The  Aircraft  Equa¬ 
tions  of  Motion  Block,  constructed  by  Kuechenmeister  [Ref.  2],  incorporates  a 
nonlinear  six  degree  of  freedom  equations  of  motion  model.  The  stability  dmva- 


tives  used  in  the  model  are  constants  either  calculated  or  estimated  for  a  tK>minal 
cruise  flight  condition  of  73  feet  per  second  airspeed  at  sea  level  altitude. 

The  sensor  motkls  developed  in  [Ref.  2]  were  obtained  from  manufacturing 
data  used  in  conjunction  with  the  outfitting  of  an  actual  NPS  UAV,  Archytas,  with 
accelerometers,  rate  gyros,  and  inclinometers.  Errors  introduced  into  the  sensors 
were  obtained  from  the  manufacturer’s  data  for  the  rate  gyros  and  accelerometers 
and  from  [Ref.  1]  for  the  inclinometers  and  payload  sensors .  These  errors  were 
specified  for  the  modular  integrated  avionics  group  (MIAG)  architecture  in  the 
report  and  represent  the  most  sensitive  avionics  and  navigation  system  currently 
available  for  an  unmanned  aerial  vehicle. 

The  Equations  of  Motion  Block  solves  the  nonlinear  equations  for  nine  states; 
three  body-fixed  translational  velocities,  three  body-fixed  rotational  velocities,  and 
three  Euler  angles  used  for  rotation  between  a  body-fixed  and  an  inertial  (earth- 
fixed)  coordinate  system.  The  development  of  these  equations  of  motion  is  well 
covered  in  [Ref.  2].  The  UAVs  true  path  over  the  ground  is  determined  by  resolv¬ 
ing  body  acceleration  into  inertial  (earth-fixed)  coordinates  via  a  coordinate  trans¬ 
formation  and  integrating  the  acceleration  twice  to  obtain  displacement  in 
Cartesian  coordinates.  The  Equations  of  Motion  Block,  since  it  calculates  both  the 
aircraft  states  and  inertial  position  at  each  point  in  time,  is  the  source  of  uncor¬ 
rupted  navigational  and  attitudinal  data  used  as  the  'Truth  Model"  in  the  Monte 
Carlo  simulation  and  analysis,  which  is  discussed  in  Chapter  IV. 

The  Sensor  Block  models  accelerometers,  rate  gyros,  and  inclinometers.  These 
"measure"  body-fixed  accelerations,  body  angular  rates,  and  Euler  Angles  respec¬ 
tively.  The  models,  covered  in  [Ref.  2],  incorporate  both  biases  and  sensor  noise 


floor  errors.  Each  sensor  is  modeled  as  a  flrst-order  filter  with  a  manufacturer-sup¬ 
plied  cutoff  frequency.  The  accelerometer  and  rate  gyro  anti-aliasing  Alters  were 
modeled  as  first  order  low-pass  Alters. 

The  Sensor  Block  ouq)uts  "measured”  body-Axed  accelerations,  body  angular 
rates,  and  Euler  Angles.  These  measured  parameters  are  passed  to  the  Navigation 
Block. 

B.  NAVIGATION  SYSTEM 

The  Navigation  Block  incorporated  two  different  but  complementary  naviga¬ 
tion  systems:  an  inertial  navigation  and  Global  Positioning  System  (GPS).  While 
GPS  is  setting  new  standards  for  navigational  position  accuracy,  its  slow  update 
rate  (once  per  second  at  best)  makes  it  less  than  perfectly  suited  for  rapid  position 
updates  for  a  fast-moving  aircraft.  On  the  other  hand,  an  inertial  navigation  sys¬ 
tem,  which  measures  body-Axed  accelerations,  converts  them  via  coordinate  trans¬ 
formation  to  inertial  accelerations  and  integrates  them  to  determine  velocity  and 
position,  is  good  in  the  short  term  but  tends  to  introduce  position  errors  over  time. 
These  errors  are  a  result  of  small  biases  in  acceleration  sensors.  One  answer  to  the 
problem  is  to  re-initialize  inertial  position  as  often  as  possible.  A  dynamic  vehicle 
like  a  UAV  requires  position  updates  many  times  a  second,  a  capability  inherent  in 
inertial  navigation  systems.  The  solution  to  providing  quick  position  updates  and 
small  steady  state  errors  is  to  design  a  complementary  Alter  which  incorporates  the 
best  features  of  both  GPS  and  INS. 


The  navigation  system  designed  in  [Ref.  3]  incorporates  a  strapdown  inertial 
navigation  system  consisting  of  accelerometers  measuring  acceleration  along  the 
aircraft  x  (positive  forward),  y  (positive  toward  the  right  wing-tip),  and  z  (positive 
down)  axes  and  gyros  measuring  angular  rates  p  (about  the  aircraft  x  axis),  q  (about 
y),  and  r  (about  z).  Inclinometers  in  the  sensor  module  are  used  to  measure  aircraft 
roll  (phi),  pitch  (theta),  and  yaw  (psi)  angles.  The  Euler  angles  phi,  theta,  and  psi 
are  used  to  transform  the  accelerations  sensed  in  the  body  coordinate  system  to  an 
earth-fixed  inertial  coordinate  system.  The  complete  navigation  system  [Ref.  3] 
incorporated  a  Kalman  filter  to  fuse  the  low-frequency  GPS  updates  with  the  higher 
frequency  accelerometer/rate  gyro  acceleration  data.  To  imixove  accuracy,  the  sys¬ 
tem  included  a  differential  GPS  model.  This  incorporates  a  ground-based  system 
placed  on  a  surveyed  site  which  broadcasts  its  measured  GPS  error  ro  the  local  area 
[Ref.  3]  and  is  used  to  improve  the  accuracy  of  the  onboard  GPS  system. 


C.  GUroANCE 

The  guidance  block  in  Figure  2  was  developed  as  a  part  of  this  work  and  took 
the  filtered  X  and  Y  inertial  position  data  trom  the  navigation  block  and  compared 
it  with  the  desired  inertial  position.  More  details  are  included  in  the  next  chapter. 


D.  CONTROLLER 

The  controller,  developed  by  the  author,  was  implemented  to  control  four  vari¬ 
ables  in  steady  state-airspeed,  attitude,  ground  track  heading,  and  rudder  position. 
It  is  discussed  in  full  in  the  following  chapter. 
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E.  PAYLOAD  MODULE 


The  payload  module  in  Figure  2  simulated  a  low-light  television  camera  or 
FLIR  used  to  target  an  enemy  position.  Because  there  was  no  feedback  from  the 
payload  module  to  the  rest  of  the  model  and  consequently  no  jn-opagation  of  pay- 
load  sensor  errors,  the  payload  was  modeled  as  a  sensor  with  simple  lookdown  and 
heading  errors  as  given  for  the  MIAG  payload  in  the  UAV  Error  Analysis  Report 
[Ref.  1].  There  is  currently  no  servo  model  included  in  the  payload  module. 
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ra.  CONTROLLER  DESIGN,  GUIDANCE,  AND 

PAYLOAD  MODEL 


The  controller  was  designed  using  Linear  Quadratic  Regulator  techniques.  The 
goal  was  to  construct  a  closed-loop  model  which  jrielded  satisfactory  performance 
while  quickly  damping  out  transient  responses  so  that  steady-state  errors  could  be 
easily  evaluated  by  Monte  Carlo  analysis  of  a  simulated  flight  The  notation 
adofHed  here  means  to  be  suggestive.  Upper  case  letters  are  used  to  denote  the  full- 
scale  values  of  aircraft  variables  or  matrices,  while  lower  case  letters  indicate  per¬ 
turbations  around  the  trim  values  of  the  full-scale  variables.  A  subscripted  "B” 
indicates  a  variable  exi»essed  in  body  coordinates  while  a  subscripted  "U"  indi¬ 
cates  one  exfuessed  in  inertial  coordinates.  A  coordinate  transformation  from  iner¬ 
tial  to  body  is  indicated  by 


and  a  body  to  inertial  transformation  is  shown  by  exchanging  the  superscript  and 
subscript 


A.  PLANT  DESCRIPTION 

Since  the  controller  design  required  a  linear  plant  the  equations  of  motion 
model  [Ref.  2]  was  linearized  at  a  cruise  condition  of  73.3  feet  per  second  true  air¬ 
speed  at  a  sea  level  standard  day.  The  linearized  model  has  the  following  form: 
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X  =  Ax  +  Bu 
y  =  Cx  +  Du 


(EQl) 


where  A  and  B  are  given  in  Appendix  B,  C  is  a  9x9  identity  matrix,  and  D  is  a  4x9 
zero  matrix.  The  state  vector  x  consisted  of  the  following  states: 

X  =  ja^  u  V  w  p  q  r  <1)  0  Xjfj'.  (eq2) 

where  the  first  four  states  were  the  elevator,  rudder,  aileron,  and  throttle  actuator 
portions.  The  states  u,  v,  and  w  were  linear  velocities  about  the  aircraft  x  (positive 
forward),  y  (positive  out  the  tight  wing)  and  z  (positive  down)  axes  respectively. 
The  three  states  p,  q,  and  r  were  rotational  velocities  about  the  aircraft  x,  y,  and  z 
axes  which  used  the  right-hand  rule  to  determine  their  orientation.  The  final  three 
states,  6,  and  tjr,  were  Euler  angles  used  to  determine  the  aircraft’s  orientation 
with  respect  to  an  earth-fixed  coordinate  system  [Ref.  2].  The  control  inimt  vector 
u  is  given  by: 


u  = 


(EQ3) 


where  the  first  dtree  elements  in  u  represented  deflections  in  elevator,  rudder,  and 
aileron  control  surfaces.  The  fourth  element,  S^,  was  measured  throttle  input  as  a 
percentage  of  total  thrust  available.  The  output  vector  y  was  the  set  of  the  mea¬ 
sured  outputs  obtained  from  the  sensor  suite.  The  onboard  sensors  could  not 
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directly  measure  all  of  the  states  as  formed  in  the  original  model;  therefore,  the 
original  set  of  nine  states  was  replaced  by  nine  measurable  outputs  y  : 


y  = 


a  a  a  a  v 
e  r  a  t  t 


“yijPqf't'eVg, 


(EQ4) 


The  vector  y  was  obtained  by  premultiplying  the  original  set  of  states  x  by  a  matrix 
T: 

y  =  Tx  (EQ5) 

Since  the  transformation  matrix  T  is  composed  of  constant  elements  taken  about 
the  cruise  trim  condition,  we  can  also  write 
y  =  Tx. 

Substituting  into  Eq.  1,  we  obtain 


T  V  =  AT  ^y  +  Bu.  (EQ6) 

In  state  space  form,  Eq.  6  becomes 

y  =  TAT  ^y  +  TBu.  (EQ7) 

The  matrix  T  is  called  a  similarity  transformation.  Eigenvalues  for  the  trans¬ 
formed  plant  matrix  T'^  AT  (which  are  the  same  as  the  eigenvalues  of  the  original 
plant  A)  are  (resented  in  Table  1 .  The  linearized  model  was  stable  in  all  modes 
with  the  exception  of  a  slightly  unstable  spiral  mode. 
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The  ouqnits  differed  from  the  original  states  in  the  following  manner; 

1.  Forward  velocity  u  was  replaced  by  true  airspeed  v^. 

2.  Lateral  and  vertical  velocities  v  and  w  were  replaced  by  lateral  and  vertical 
accelerations  ny  and  n2. 

3.  Heading  \|/  was  replaced  by  ground  track  heading 


TABLE  1:  BLUEBIRD  OPEN  LOOP  EIGENVALUES 


Mode 

Eigenvalue 

Damping 

Short  Period 

-3.9173+/-3.4918i 

74.6% 

Phugoid 

-.0057+/-.5016i 

1.14% 

Spiral 

.0792 

0 

Dutch  Roll 

-.5322+/-3.5719i 

14.7% 

RoU 

-5.5165 

100% 

By  including  the  actuator  states  in  the  linear  model,  it  became  possible  to  use 
measurable  outputs  which  were  independent  of  actuator  inputs.  The  transformation 
matrix  is  (tefined  later  in  this  chapter.  Since  there  were  the  same  number  of  mea¬ 
sured  ouQiuts  and  (Miginal  states,  it  was  possible  to  design  a  state-feedback  control¬ 
ler  for  this  set  of  ou^ts. 

The  first  transformation  was  effected  because  the  pitot-static  system  aboard 


Bluebird  did  not  measure  velocity  along  the  aircraft  x-axis  but  rather  measured  the 
velocity  in  wind  coordinates.  This  is  resolved  in  the  cruise  condition  that  the  con¬ 
troller  was  designed  around  by: 


(EQ8) 


COS0Q  sin0Q 

U 

L.  .« 

w 

where  @o  is  the  trim  aircraft  pitch  attitude  of  .0912  radians  (api»x>ximately  five 
degrees).  Sideslip  (v)  and  body  vertical  velocity  (w)  can  be  obtained  practically 
only  by  integrating  body  y  and  z  accelerations  ny  and  n^,  which  are  available 
directly  from  the  accelerometers  of  the  strapdown  inertial  navigation  system. 
Therefore,  v  and  w  were  replaced  with  ny  and  n^.  This  was  accomplished  by 
using  the  second  and  third  rows  of  the  original  plant  (A)  matrix,  since  in  the  linear 
form,  state  derivatives  can  be  formed  by  multiplying  the  second  and  third  rows  of 
the  original  plant  matrix  by  the  state  vector. 


(EQ9) 


where  A(2)  and  A(3)  represent  the  second  and  third  rows  of  the  original  plant 
matrix.  The  fourth  and  final  transformation  was  from  heading  (y)  to  ground  track 
heading  (tjfgt).  The  change  to  ground  track  heading,  one  of  four  variables  con¬ 
trolled  in  steady  state,  was  necessary  to  keep  the  UAV  ground  track  always  pointed 
in  the  direction  of  the  desired  waypoint  Without  this  transformation,  the  guidance 
routiiw  would  have  continually  commanded  aircraft  leading  toward  the  targeted 
waypoint  Any  crosswind  would  cause  a  cross-track  error,  with  the  aircraft  follow¬ 
ing  a  banana-shaped  route  to  the  taigeted  waypoint  With  the  transformation  to 
ground  track  heading  and  the  controlling  of  rudder  position  to  zero  in  steady  state. 
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the  aircraft  naturally  turns  into  the  wind  and  flies  a  direct  path  to  the  next  waypoint. 
Fortunately,  ground  track  heading  is  a  combination  of  states  [Ref.  4]; 


B-Asina 

\j/+  - 

cosy 


(EQ  10) 


Equation  10  is  derived  in  Appendix  A.  The  newly  introduced  variables  p,  y,  and  a 
can  be  computed  using  their  steady  state  values  at  the  design  cruise  condition. 
Assuming  a  small  sideslip  angle  in  cruise,  ^  (sideslip)  can  be  approximated  by 


(EQ  11) 


Since  the  vehicle  is  trimmed  in  a  straight  and  level  cruise  condition,  y  (the  dif¬ 
ference  between  angle  of  attack,  cx,  and  inertial  flight  path  angle,  6)  is  negligible. 
For  the  design  cruise  condition,  with  the  aircraft  in  Ig,  non-climbing  flight,  the  lift 
L  generated  by  the  aircraft  must  be  equal  to  the  weight  of  the  vehicle.  Using  vehi¬ 
cle  parameters  [Ref.  2]  and  a  relation  for  steady-state  lift  from  Lan  and  Roskam 
[Ref.  5], 


2 

L  =  0.5pvV.-sCj  a  =  Weight,  (eq  i2) 

where  p  is  air  density,  s  is  wing  area,  0^^  is  the  aircraft  lift-curve  slope,  and  a  is 
angle  of  attack.  Using  the  Bluebird  parameters  from  Table  1  for  the  design  cruise 
condition  at  standard  temperature  and  pressure. 
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a  =  0.0912. 


(EQ  13) 


Thus,  Eq.  7  simplifies  to 


Now,  the  measured  ouq)uts  can  be  derived  from  the  original  state  vector  using  the 
9x9  transformation  matrix  T: 


u 

< 

"y 

P 

q 

r 

e 

y. 

cos  (Gq) 
A  (2) 
A  (3) 

0 

0 

0 

0 

0 

0 


0  sin  (Gq)  0  0  0  0 


OOi 


0 

0 

0 

0 

0 

2 

73.3 


0 

0 

0 

0 

0 


1  0  0 
0  1  0 
00  1 
00  0 
000 


0 

0 

0 

1 

0 


— -  0  0  0  0  -sin  (0.09)  0  1 


) 

j- 

u 

V 

) 

w 

) 

p 

) 

q 

) 

r 

3 

e 

y 

(EQ15) 


B.  DESIGN  REQUIREMENTS 

The  controller  was  required  to  make  all  closed  loop  eigenvalues  of  the  Bluebird 
model  stable  and  meet  a  minimum  requirement  of  70%  damping  and  a  maximum 
eigenvalue  frequency  of  12  radians  per  second.  The  damping  requirement  was 
mandated  to  prevent  overshoot  in  system  response  to  steady-state  commands  while 
the  maximum  eigenvalue  requirement  was  present  in  consideration  of  actuator 
bandwidth  limits.  Because  actuators  were  limited  in  frequency  response  to 
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12  rad/s,  any  system  eigenvalues  faster  than  this  would  have  resulted  in  modes  too 
fast  for  the  control  surfaces  and/or  throttle  to  deal  with  [Ref.  6]. 


C.  CONTROL  STRATEGY 

Since  only  four  controls  (elevator,  aileron,  rudder,  and  throttle)  were  available, 
only  four  variables  could  be  controlled  in  steady  state.  These  controlled  states  were 
airspeed,  pitch  angle,  ground  track  heading,  and  rudder  position.  Two  of  these 
states  were  in  the  aircraft  x-y,  or  lateral,  plane  (ground  track  heading  and  rudder 
position),  and  two  are  in  the  x-z,  or  longitudinal,  plane  (airspeed  and  theta).  This 
selection  of  two  states  in  each  plane  was  made  because  two  of  the  controls,  aileron 
and  rudder,  are  largely  lateral  surfaces,  and  the  other  two,  elevator  and  throttle,  are 
largely  longitudinal.  In  the  design  cruise  condition,  there  is  nearly  zero  coupling 
between  the  two  planes.  This  will  be  demonstrated  by  the  controller  feedback  gain 
matrix. 

In  order  to  keep  the  nose  of  the  aircraft  turned  into  the  wind  when  flying  its 
command  ground  track  heading,  rudder  position  was  selected  as  one  of  the  con¬ 
trolled  states.  Because  its  augmented  plant  is  laterally  stable.  Bluebird  is  guaran¬ 
teed  to  have  a  finite  steady-state  sideslip  angle.  The  fact  that  Bluebird  (and  for  that 
matter  the  vast  majority  of  aircraft)  is  symmetrical  in  the  body  x-z  plane  guarantees 
that  this  sideslip  angle  will  be  zero.  In  order  to  maintain  this  symmetrical  condi¬ 
tion,  the  rudder  must  be  commanded  to  zero  in  steady  state  [Ref.  4].  The  other  con- 
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trolled  lateral  state,  ground  track  heading  (Eq.  10),  ensures  that  the  aircraft  flies  a 
heading  in  the  inertial  X-Y  plane  toward  the  next  selected  waypoint. 


D.  SYNTHESIS  MODEL 

The  controller  synthesis  model,  shown  in  Figure  3,  was  formed  by  appending 
the  derivative  and  error  outputs  for  the  variables  controlled  in  steady  state  to  the 
transformed  Bluebird  linear  model.  It  served  as  an  interface  between  the  designer 
and  control  algorithm  formed  using  the  linearized  Bluebird  plant  with  die  trans¬ 
formed  state  vector  previously  described.  To  the  synthesis  model  were  added  the 
command  inputs  (rudder  position,  airspeed,  ground  track  heading,  and  attitude). 
Command  errors  were  formed  by  subtracting  aircraft  states  from  commanded  states 
and  integrators  were  placed  on  these  errors  to  drive  diem  to  zero  in  steady  state 
[Ref.  4].  The  creation  of  the  linear  synthesis  model  makes  possible  the  iterative 
approach  used  in  forming  the  aircraft  controller. 


E.  LINEAR  QUADRATIC  REGULATOR  DESCRIPTION 
Consider  the  following  linear  system: 


T|  =  AT|-t-Bjr-t-B2U 
z  =  CjT|  +  DjU 


(EQ16) 


where  q  is  composed  of  the  measured  ouqiuts,  rudder  position,  airspeed,  attitude, 
and  ground  track  heading  command  error  states: 
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Convnand 


Figure  3.  Bluebird  Controller  Synthe^s 
Model 


20 


Ti  =  [y  5e^  6e^  5e^  (eq  i7) 

Equation  16  is  a  state  space  re{Mnesentation  of  Figure  3.  We  assume  that  (A.Ba)  is 
stabilizable,  (C},A)  is  detectable,  and  D  is  full  column  rank.  The  m^x  B]  is  the 
command  input  matrix,  B2  is  the  control  surface  input  matrix,  and  r  is  the  set  of  ref¬ 
erence  commands  that  are  desired  to  be  tracked.  Now,  we  define  a  cost  function,: 

J  =  U‘“{CQC-»-D'RD}dt.  (EQI8) 

wh^  dis^onal  matrices  Q  and  R,  used  to  weight  the  ouq)uts  and  inputs,  ate  posi¬ 
tive  semidefinite  and  positive  definite  respectively.  Hie  Linear  Quadratic  Regula¬ 
tor  (LQR)  iMoblem  is  to  find  the  state  feedback  controller  K,  wher>; 

K=[kpKj].  (EQ19) 

such  that  the  feedback  system  shown  in  Hgure  4  is  stable  and  J  is  minimized.  It 
can  be  shown  that 


K  =  R“^B2'P,  (EQ20) 

where  P  is  a  stabilizing  solution  of  the  Algebraic  Ricatti  Equation  [Ref.  10]: 

AP  +  PA'  -  PB2R”^B' P  -h  C'QC  =  0  (EQ  21) 
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Figure  4.  Feedirack  System 


F.  CONTROLLER  DESIGN 

As  seen  in  Figure  3,  the  controller  synthesis  model  placed  an  integrator  on 
errors  between  measured  and  desired  airspeed,  theta,  ground  track  heading,  and 
rudder  position.  If  a  steady-state  error  exists  in  any  of  these  states,  then  the  integral 
of  that  error  creates  an  actuating  signal  to  drive  that  error  to  zero.  Since  four  new 
error  states  are  created  by  the  addition  of  the  integrators,  the  construction  of  tl^ 
LQR  controller  creates  the  proper  feedback  gains  for  the  error  states  as  well  as 
those  of  the  transformed  Bluebird  model  with  actuators  [Ref.  4].  To  create  a 
smoother  transient  response,  the  derivatives  of  velocity,  theta,  and  ground  track 
heading  were  also  made  outputs,  adding  them  to  the  cost  function  to  be  described 
shortly.  The  (terivative  terms  add  a  useful  tool  to  the  controller  design  process. 
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Adding  weight  to  the  state  derivative  ou^uts  (via  the  Q  matrix)  increases  system 
damping,  reducing  or  eliminating,  when  damping  exceeds  .707,  overshoot  in  sys¬ 
tem  response  [Ref.  7].  In  the  design  of  the  controller,  it  is  important  to  ensure  that 
all  eigenvalues  of  the  closed  loop  plant,  A-B2K»  are  slower  (of  lower  frequency) 
than  the  bandwidth  of  the  control  loops  (the  frequency  responses  between  the  ele¬ 
vator,  rudder,  aileron,  and  throttle  inputs  and  the  outputs  of  the  feedback  gain 
matrix  of  the  closed  loop  system)  [Ref.  6].  Also,  since  anticipated  system  com¬ 
mands  in  velocity,  theta,  and  ground  track  heading  are  under  1  rad/s,  command 
loop  bandwidths  ( the  frequency  responses  between  command  inputs  and  the  like 
state  ouqnits  of  the  closed  loop  plant)  in  the  area  of  1  radian  per  second  were  con¬ 
sidered  acceptable. 

The  Matlab  file  used  to  create  the  LQR  Controller,  EVAL.M,  is  given  in  Appen¬ 
dix  A.  To  achieve  the  specified  design  requirements,  two  tools  available  to  the  con¬ 
trol  designer  are  the  Q  and  R  matrices.  By  using  the  diagonal  elements  in  these 
matrices  to  weight  the  synthesis  model  inputs  and  outputs,  the  designer  can  tune  the 
weights  iteratively  to  obtain  desired  control  and  command  loop  bandwidths,  eigen¬ 
values,  and  system  damping.  To  begin  design  of  the  controller,  Q  and  R  are  set  as 
identity  matrices  of  size  7x7  and  4x4,  corresponding  to  the  number  of  system  out¬ 
puts  and  control  inputs  respectively.  Next,  the  maximum  (highest  frequency) 
eigenvalue  of  the  closed  loop  linearized  system  and  the  minimum  damping  associ¬ 
ated  with  the  closed  loop  eigenvalues  are  recorded.  After  this,  the  frequency 
response  of  the  command  and  control  loops  are  observed  and  the  associated  band- 
widths  are  observed.  Finally,  the  time  constants  (the  time  in  seconds  for  a  state  out- 
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put  to  reach  63%  of  a  commanded  step  input  value)  for  velocity,  theta,  and  ground 
track  headings  are  recorded  as  a  measure  of  system  performance. 

Analysis  of  the  controller  begins  with  the  control  loop  bandwidths:  the  closed- 
loop  frequency  responses  between  the  control  (elevator,  rudder,  aileron,  and  throt¬ 
tle)  inputs  and  the  feedback  gain  matrix  ouqnit.  These  bandwidths  are  obtained  by 
"cutting"  the  closed-loop  system  at  the  input  and  treating  the  output  of  the  feedback 
gain  matrix  as  the  system  output  (Figure  4).  This  makes  our  system,  for  analysis 
purposes: 

X  =  ( A  -  B2K)  X  +  B2U 

U  =  -KX  (EQ  22) 

If  any  closed  loop  eigenvalue  exceeds  the  maximum  actuator  bandwidth  (in  the 
case  of  Bluebird,  12  radians  per  second),  then  the  control  loop  with  the  maximum 
bandwidth,  which  is  generally  responsible  for  the  maximum  eigenvalue,  is  penal¬ 
ized  by  increasing  its  respective  diagonal  element  in  the  R  matrix.  For  example,  if 
rudder,  the  second  input,  had  the  widest  bandwidth,  the  the  (2,2)  element  in  the  R 
matrix  is  increased.  Unfortunately,  there  is  no  formula  for  determining  how  much 
to  penalize  a  control  input.  It  is  an  iterative  process  which  is  repeated  until  no 
eigenvalues  exceed  the  maximum  allowable  bandwidth. 

The  command  loop  bandwidths  are  made  larger  by  increasing  the  respective 
diagonal  elements  in  the  Q  matrix.  Again,  this  is  an  iterative  process  in  which 
adjustments  are  made,  the  new  command  bandwidths  are  observed,  and  then  new 
adjustments  are  made  in  Q.  The  command  loop  bandwidths  are  observed  by  look- 
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ing  at  the  response  from  the  command  inputs  (rudder,  velocity,  pitch  angle,  and 
ground  track  heading)  to  the  outputs  of  the  commanded  states  (Figure  4) 


x=  (A-B2K)x  +  B,r 

y2  =  CtX 


(EQ23; 


where  Ci  is  the  output  matrix  for  the  four  commanded  states  only. 

The  final  step  is  the  analysis  of  the  response  of  the  linear  model  to  step  inputs. 
Any  overshoot  can  be  reduced  by  increasing  weight  on  the  respective  state  deriva¬ 
tive  output,  which  controls  damping  [Ref.  6].  Obviously,  the  addition  of  too  large  a 
weight  on  the  derivatives  will  slow  the  system  response  to  command  and  control 
inputs. 


G.  BLUEBIRD  CONTROLLER  LINEAR  RESULTS 

The  design  iterations  for  the  LQR  controller  are  given  in  Table  2.  The  table 
presents  the  diagonal  elements  of  the  Q  and  R  matrices,  maximum  closed  loop 
eigenvalues,  and  control  and  command  loop  bandwidths. 


TABLE  2:  BLUEBIRD  CONTROLLER  SYNTHESIS 


Q 

R 

Max 

E- 

Control  Loop  Bandwidths 

Cinnmand  Loop 
Bandwidths 

val 

Bev 

Rud 

Ail 

Thr 

Vel 

'For 

e 

[....] 

[i . . 

-86 

30 

2 

10 

10 

1 

.2 

.8 

[.  W  ,] 

.100  1  1  1 

-44 

.10 

1.2 

10 

10 

1 

.2 

.8 

[....] 

[lOO  1  1  10^ 

-11 

10 

1.2 

10 

1 

.8 

.2 

.8 

[l  1  100  l] 

[lOO  1  1  10^ 

-11 
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The  feedback  gain  matrix  is  as  in  Eq.  19,  where  for  the  last  iteration, 


0.99  0  0  -0.00  0.08  0  0  0  -0.53  0  0  -3.3  0 

0  0.13  0.40  0  0  0.06  0  -0.07  0  3.5  -2.4  0  -1.8 

0  0.13  0.24  0  0  -0.02  0  0.1  0  -1.2  0.92  0  0.45 

-0.01  0  0  0.49  0.13  0  0  0  -0.18  0  0  -3.8  0 


(EQ24) 


and 


-0.0992  0  0  -0.1273 

0  -0.9923  0.1236  0 

0  0.1236  0.9923  -0.0001 

_0.0127  0  0  0.9919_ 


(EQ25) 


The  columns  of  each  element  of  K  are  the  feedback  gains  for  the  control  actua¬ 
tor  positions  (columns  1-4),  the  measured  ouq>uts  (columns  5-13)  and  the  com¬ 
mand  error  integrators  (columns  14-17).  The  rows  of  K  correspond  to  the  elevator, 
rudder,  aileron,  and  throttle  commands  respectively.  As  expected  for  the  design 
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cruise  coiKlition,  the  longitudinal  (elevator  and  throttle)  and  lateral  (rudder  and 
aileron)  planes  are  largely  uncoupled.  The  resulting  controller  has  the  following 
state  space  reinresentation: 

Xc  =  r-y2 

(EQ  26) 

u  =  KjX^  +  Kpy 

where  x^.  is  the  set  of  four  states  formed  from  the  command  errors. 

The  airspeed-pitch  angle  longitudinal  combination  provides  a  great  degree  of 
flexibility.  Since 

0  =  a-l-Y  (EQ27) 

and  in  the  design  condition, 

Y  =  0,  (EQ28) 

we  arrive  at: 

0  =  a.  (EQ29) 

Although  angle  of  attack  sensors  are  widely  used  and  one  is  available  on  Blue¬ 
bird,  they  are  typically  noisy  and  not  well  suited  for  feedback  [Ref.  6].  In  the  zero 
flight  path  angle  cruise  condition,  6  serves  as  a  good  approximation  for  angle  of 
attack.  From  Eq.  7,  which  is  good  for  Ig  climbing  as  well  as  straight  and  level 
flight,  for  a  constant  airspeed,  angle  of  attack  will  be  constant.  Any  increase  in  6 
will  be  reflected  in  a  like  increase  in  y.  If  we  assume  a  wings-level  condition,  then 


(EQ30) 


V^siny  =  z 

~  where  z  is  inertial  altitude.  From  this  relation,  a  simple  climb  rate  controller 
could  be  devel(^)ed.  To  preserve  simplicity  and  reduce  computation  time,  this  was 
not  imi^emented. 

The  control  loop  response  plots  are  shown  in  Figures  5-8.  All  control  loop 
bandwidths  were  kept  below  the  12  rad/s  limit.  The  line  which  runs  horizontally 
near  0  dB  to  the  respective  bandwidth  is  the  main  control  loop  response  while  the 
<^her  three  lines  in  each  plot  are  the  responses  of  each  of  the  three  other  controls  to 
die  input  control. 

The  command  loop  response  plots  are  shown  in  Figures  9-11.  The  three  com¬ 
mand  loop  bandwidths  were  placed  near  die  target  of  1  rad/s.  Linear  model 
responses  to  step  inputs  in  velocity,  ground  track  heading,  and  theta  commands  are 
shown  in  Figures  12-14.  Time  constants  (times  to  reach  63%  of  steady  state  value) 
for  these  responses  are  9,  4,  and  8  seconds  and  minimal  overshoot  is  present  for 
each  command  response. 


H.  IMPLEMENTATION  ON  THE  NONLINEAR  MODEL 

The  complexity  of  the  complete  Bluebird  model  with  sensor  and  GPS/INS  nav¬ 
igation  models  precluded  linearization  with  the  LQR  controller  in  the  loop  to  deter¬ 
mine  system  eigenvalues  and  prove  stability.  While  the  extensive  system 
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EI«v«itor  Comnn«nol*d  to  EI«v«tor 


Figure  5.  Elevator  Commanded  to  Elevator 


Ruddor  ComrTMiirtdod  to  Ruddor 


Figure  6.  Rudder  Commanded  to  Rudder 
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AII«ron  CommatrtclAcd  to  Alloron 


Figure  7.  Aileron  Conunanded  to 
Aileron 


XMrottI*  Oommandcict  to  Xt^rottl* 


Figure  8.  Throttle  Comnuuided  to 
Ilirottie 
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\/9loclty  Corrtmetnciscl  to  N/elocIty 


Orourtd  Xracic  CommandACl  to  Oround  Xi-acK 


Figure  11.  Ground  'nuck  Conunanded  to  Ground  IVack 


RESPONSE  TO  A  STEP  Vdloelty  OOIS4MAND 


Figure  12.  Response  to  a  Step  Velocity  Command 
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RESRONSE  TO  A  SXER  ATTITUDE  COMMAND 


Figure  13.  Response  to  a  Step  Attitude 
Conunand 


RESRONSE  TO  A  STER  REAOINO  COMMAND 


Figure  14.  Response  to  a  Step  Heading 
Conunand 


33 


dynamics,  including  sensor  models  and  the  navigation  Kalman  niter,  may  be  de¬ 
stabilizing,  the  best  available  indication  of  plant  stability  is  that  the  model  always 
returned  to  a  stable  cruise  condition  even  when  simulation  was  begun  substantially 
off  of  the  design  cruise  condition. 

The  states  used  in  feedback  were  "measured”  using  the  sensor  models  [Ref.  2  ] 
and  the  navigation  filter  [Ref.  3].  The  airspeed  sensor  used  the  aircraft  forward 
velocity,  resolved  in  wind  coordinates,  with  zero-mean  white  noise  added  to  give 
an  RMS  error  of  .5  ft/s.  The  accelerometer  models  used  to  measure  acceleration 
incorporated  both  gaussian  zero-mean  errors  reflecting  the  modeled  sensor  noise 
floors  as  well  as  cross-axis  errors  equal  to  0.5%  of  the  off-axis  acceleration.  The 
rate  gyro  models  also  incorporated  a  noise  floor  with  a  0.5%  cross-axis  error.  The 
inclinometers  incorporate  a  0. 1  degree  RMS  error  with  a  0.5%  cross-axis  error.  All 
sensors  were  modeled  as  first-order  filters  with  pass  bands  as  specified  by  the  man¬ 
ufacturer  [Ref.  2].  The  navigation  filter  was  used  to  calculate  inertial  velocity  as 
well  as  inertial  position,  which  was  fed  to  the  guidance  module  for  calculation  of 
ground  track  heading  command.  Velocity  and  theta  commands  were  open  inputs 
which  could  be  changed  as  the  simulation  was  running. 

The  controller  was  im{demented  on  the  nonlinear  model  using  Delta  implemen¬ 
tation  [Ref.  8].  This  involved  taking  the  feedback  integrator  that  was  placed  on  the 
command  errors  (to  drive  them  to  zero  in  steady  state)  and  moving  it  to  the  input  of 
the  nonlinear  plant.  To  maintain  consistency,  state  derivatives,  rather  than  the 
states  themselves,  were  used  in  feedback;  the  state  derivatives,  multiplied  by  the 
constant  feedback  gains  were  now  integrated  at  the  plant  input.  Recall  that  the  con¬ 
troller  previously  obtained  had  the  form: 
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r  ■■  ■ 

ii':  • 

I 


\  =  *'~y2 

(EQ30) 

u  =  KjX^  +  Kpy 


By  moving  the  integrator  on  the  command  error  to  the  plant  input,  the  linear 
controUo'  becomes: 


(EQ31) 


Figure  IS  shows  the  plant  with  the  Delta-imj^emented  controller. 


Figure  IS.  Contitrtler  With  Delta  Implementation 


Further  ramifications  of  Delta  implementation  are  discussed  in  [Ref.  8].  For  the 
puiposes  of  this  controller  design,  the  pl£K;ement  of  an  integrator  at  the  plant  input 
provided  an  easy  way  to  set  the  control  surface  initial  conditions  to  their  trimmed 
cruise  values.  This  was  important  in  reducing  startup  transients  when  the  Simula- 
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tion  was  run;  limits  on  memory  and  time  required  that  as  much  of  each  run  as  pos¬ 
sible  accurately  reflected  a  model  in  stable  flight 


I.  GUIDANCE 

The  guidance  mechanism  stored  three  waypoints  which  were  selected  to  navi¬ 
gate  the  model  around  a  fixed  target.  The  waypoint  guidance  used  the  error 
between  position  estimated  in  the  navigation  filter  and  the  current  waypoint  to 
command  a  ground-track  heading  to  the  waypoint  This  heading  command  was 
found  by  first  taking  the  differences  AX  and  AY  between  the  estimated  and  the 
desired  inertial  position  at  each  time  step: 

AX  =  X^p  - 

(EQ32) 

AY  =  Y^p  -  Y^g^j. 

with  X^  and  Y^p  denoting  the  waypoint  coordinates  and  and  YgsT^ing 

the  position  estimates  from  the  navigation  block.  Next  the  ground  track  heading 
was  computed  by  solving  the  arctangent  of  the  error  components: 

AY 

¥gt  =  atan— .  (EQ33) 

When  Bluebird  came  within  two  seconds  of  a  selected  waypoint  logic  within  the 
guidance  block  commanded  the  next  waypoint  Because  of  simulation  time  con¬ 
straints,  only  three  waypoints  were  stored.  Appendix  A  contains  the  Matlab  code 
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for  the  guidance  routine,  PSICOM2.M.  The  Simulink  guidance  block  is  shown  in 
Hguie  16. 


Pointai 


Figure  16.  Goidance  Block 


J.  PAYLOAD 

The  simulated  payload  was  modeled  as  an  open-loop  sensor  with  azimuth  and 
lookdown  errors  of  0..285  degrees  at  one  standard  deviation  (la).  Because  there 
was  no  feedback  from  the  payload  sensors  to  the  aircraft,  the  payload  could  be  con¬ 
structed  as  a  simple  model  without  loss  of  fidelity  for  the  entire  simulation. 

The  camera  azimuth  and  lookdown  angles  were  found  in  three  steps.  First,  the 
differences  between  target  position  and  estimated  aircraft  position  in  the  inertial 
coordinate  position  (AX,  AY,  and  AZ)  were  found.  Next,  the  inertial  differences 
were  converted  to  body  differences  via  a  coordinate  transformation  [Ref.  2]; 
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AXg 

-  4 

B* 

AXjj 

(EQ  34) 

Ax„ 

Finally,  the  body  x,  y,  and  z  differences  were  converted  to  azimuth  (e^z)  and  look- 
down  (0ld)  angles  by: 


=  atan - 

AZ  Ax 


B 


®LD  = 


Ax^  +  Ay^ 


Az 


B 


(EQ35) 


The  Payload  Block  is  shown  in  Figure  17. 
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Figure  17.  Payload  EHagram 
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IV.  MONTE  CARLO  ANALYSIS 


A.  MONTE  CARLO  OVERVIEW 

The  teim  Monte  Carlo  was  first  used  to  describe  [Mob^ilisdc  mathematical 
method  by  scientists  working  on  the  Manhattan  Project  in  the  New  Mexico  desm 
in  die  early  1940’s.  The  mediod  is  essentially  a  series  of  games  of  chance  (or  ran¬ 
dom  events)  whose  behavior  and  outcome  can  be  used  to  study  and  classify  mathe¬ 
matical  properties  of  a  random  variable  or  a  random  ivocess.  By  repeatedly 
simulating  a  system  subject  to  noise,  the  statistical  properties  of  the  ouQHits  of  die 
system  can  be  derived.  Since  these  (xoperdes  are  derived  ftom  a  finite  number  of 
trials,  they  are  only  an  approximation  of  the  real  statistical  properties  of  the  system. 
As  might  be  expected,  these  s^^iroximations  improve  with  an  increasing  number  of 
trials. 

B.  A  BRIEF  REVIEW  OF  STATISTICS 

The  purpose  of  this  section  is  not  to  give  an  exhaustive  coverage  of  pobability 
and  statistics,  a  very  large  topic.  Rather,  a  brief  overview  of  mathematical  proper¬ 
ties  used  in  the  simulation  of  Bluebird  is  presented.  For  a  more  rigorous  coverage 
of  the  topic,  see  [Ref.  9]-[Ref.  12]. 


1.  ProlMri^ty 

Given  some  random  event  A  with  a  countable  set  of  random  outcomes 
Ai,A2,...An,  there  is  associated  with  each  mitcome  Aj^  a  {vobability  pi^  between  0 
and  1: 


0  <  pk  <  1  .  (EQ  36) 

If  the  kth  outcome  almost  never  occurs,  if  it  nearly  always  occurs,  pi^l  [Ref. 
9].  Note  the  use  of  the  word  "almost".  In  a  stochastic  process,  which  nearly  every¬ 
thing  is,  there  are  no  absolutes.  A  common  notation  for  the  probability  of  a  random 
event  A^  is 


P(A,^)  =  P|^  (EQ37) 

For  every  outcome  A^  there  is  an  associated  numerical  value  x^.  The  func¬ 
tion  assignment  X|^  to  A|^  is  termed  a  random  variable.  A  random  process  is  a  set  of 
random  variables  indexed  by  time.  The  random,  or  stochastic,  process  analyzed  in 
the  simulation  is  the  identification  by  geographical  coordinates  of  a  fixed  target  by 
the  UAV.  The  expected  value  E(x)  of  a  random  variable  x  is  defined  as 

E(X)  =  ]^piXi  =  11.  (EQ38) 

i 

The  exp^ted  value  may  also  be  thought  of  as  the  statistical  mean  of  the  random 
variable. 
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The  nth  central  moment  of  x  is  defmed  as  the  expected  value  of  the  nth 
power  of  x: 

E(  (X-H)")  =  ]^p.(x.-E(x))"  (EQ39) 

i 

The  second  central  moment  of  x  is  called  the  variance  of  x  and  is  derived 
fix)m  Eq.  40, 

E((X-li)^)  =  E((X-E(X))^).  (EQ40) 

Substituting  Eq.  40  into  Eq.  42,  we  obtain 

]^p.x?  -  E  (x)  ^  =  E  (x^)  ~  E  (x)  (EQ 41) 
i 

The  square  root  of  the  variance,  a,  is  called  the  standard  deviation  of  x  and 
is  a  measure  of  the  dispersion  of  the  random  variable  from  the  expected  value. 

A  continuous,  normally  distributed  random  variable  x,  represented  graph¬ 
ically  by  the  familiar  "bell  curve",  has  the  following  probability  density  function 
(pdf): 

f(x)  =  ,  exp - ^  (EQ42) 

aj2n  2a^ 
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Because  eirors  (Hopagate  through  the  Bluebird  simulation  in  a  complex,  nonlinear 
way,  there  are  no  guarantees  as  to  the  distributions  of  the  sensor  errors  in  the  mottel. 
However,  mean  and  standard  deviation  are  values  which  are  independent  of  distri¬ 
bution  and  so  may  be  used  to  describe  errors  in  the  Bluebird  simulation  regardless 
of  their  distribution.  Another  term  used  in  this  thesis  is  the  (Hob^iUty  distribution 
function  (PDF).  The  PDF  gives  the  probability  that  a  random  selection  of  x  will  be 
less  than  a  specified  value.  If  the  PDF  is  differentiable,  then  the  pdf  is  defined  as: 

f(x)  =-^F(x)  >0  (EQ43) 

dx 

2.  Monte  Carlo  Fundamentals 

Monte  Carlo  analysis  is  centered  around  the  idea  that  integrals  with  infinite 
limits  may  be  approximated  by  finite  sums  with  a  sufficiently  large  sample  size. 
For  a  finite  sum  of  length  N  and  distribution  G,  a  sampled  estimate  of  the  mean 
value  of  the  random  variable  x,  is  given  by 


(EQ44) 

As  N-->oo,  this  distribution  approaches  that  of  E(x).  The  observed  G  is  within  one 
standard  error  of  E(x)  68.3%  of  the  time,  two  standard  errors  95.4%  of  the  time, 
and  three  standard  errors  99.7%  of  the  time.  This  property  is  known  as  the  Central 
Limit  Theorem  of  probability  and  is  substantially  satisfied  when 


f(G)  = 


exp 


N(G-li)' 


iKf 


a 


(EQ45) 


In  otter  words,  for  a  sufficiently  lai:ge  random  sample  N,  the  [xobability  density 
function  of  the  set  of  random  samples  will  closely  approximate  the  probability  den¬ 
sity  function  of  the  continuous  random  process.  [Ref.  8] 
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V.  SYSTEM  SIMULATION 


A.  THE  SCENARIO 

The  simulation  involved  a  reconnaissance,  surveillance,  and  target  acquisition 
scenario.  The  model  was  fnst  flown  toward  a  fixed  ground  target;  acquisition  was 
at  one  nautical  mile.  Next,  a  left  turn  was  performed  in  order  to  break  closure  with 
the  target  followed  by  a  right  turn  to  keep  the  UAV  in  close  proximity  to  the  target 
For  this  scenario,  the  range  to  the  target  varied  between  3000  and  6000  feet.  In  the 
data  analysis,  all  targeting  errors  dependent  on  angular  errors  have  been  normalized 
by  distaru:e  to  the  target  to  make  them  a  furu;tion  of  angle  only.  The  Bluebird 
model  was  "flown"  over  a  set  of  three  waypoints  that  corresponded  to  the  desired 
track  with  respect  to  the  target  (Figure  18).  Superimposed  upon  the  trajectory  in 
Figure  18  are  the  inertial  X-Y  position  estimates  from  the  INS/Differential  GPS 
navigation  filter.  As  the  majority  of  the  position  error  is  in  the  inertial  vertical  (Z) 
direction,  the  X-Y  estimates  follow  the  aircraft  track  very  closely.  Although  the 
simulations  were  started  near  the  design  cruise  condition,  the  first  three  seconds  of 
each  run  was  not  included  in  the  evaluation  in  order  to  remove  any  transient  system 
response  from  the  analysis. 

Each  sensor  was  driven  by  zero  mean  band-limited  white  noise  with  a  standard 
deviation  obtained  either  from  [Ref.  1]  or,  if  it  was  not  provided  there,  from 
[Ref.  2]. 
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Overhead  View  Of  Aircraft  Trajectory 


Figure  18.  Model  IVack  Ground 


B.  NUMERICAL  INTEGRATION  PARTICULARS 

The  Runge-Kutta  S  numerical  integration  method  was  used  for  the  simulation. 
Although  Simulink  also  provides  a  more  ^curate  Adams-Gear  integration  method, 
use  of  Runge-Kutta  increased  the  simulation  speed  by  at  least  300%.  An  integra¬ 
tion  step  of  .01  seconds  was  selected  in  order  to  balance  system  performance  with 
fidelity.  At  this  time  interval,  all  noise  sources  in  the  system  had  to  be  band-limited 
to  SO  Hertz  in  order  to  satisfy  the  Nyquist  sampling  rate: 


5 


1 


(EQ46) 


2T  <- 
s  f 

max 

where  is  the  sample  time  in  seconds  and  is  the  maximum  frequency  in 
Hertz.  Selection  of  a  smaller  time  step  size  would  have  led  to  longer  simulation 
time  and  a  shorter  length  of  data;  with  the  .01  second  step,  the  computer  ran  out  of 
memory  after  approximately  ISO  seconds  of  real  time,  or  12.5  hours  of  simulation 
time. 


C.  DATA  COLLECTION 

Data  collected  for  each  run  included  the  uncorrupted  set  of  states  at  .01  second 
intervals,  the  measured  state  errors  (the  difference  between  the  actual  aircraft  states 
and  the  measured  states  used  in  feedback),  the  actual  inertial  X,  Y,  and  Z  coordinate 
position,  the  estimated  position  from  the  navigation  filter,  and  the  payload  sensor 
errors.  The  error  mean  values  and  standard  deviations  were  calculated  from  the 
measured  error  data.  Next,  the  target  tracking  error  was  calculated  at  each  point  in 
time  for  the  target  using  the  algorithm  in  Appendix  A.  The  standard  deviation  and 
mean  value  for  the  error  magnitude  were  determined  across  the  75  seconds  of  each 
run. 


D.  ANALYSIS 

The  sensor  errors  are  divided  between  the  input  errors  for  each  sensor,  assumed 
in  this  case  to  be  the  manufacturer’s  advertised  reliability  of  the  sensor  in  a  steady- 


state,  standalone  sense,  and  the  measurement  errors  obtained  from  the  Monte  Carlo 
analysis  of  the  entire  system  as  flown  around  the  set  of  waypoints.  Measurement 
errors  were  higher  in  the  simulation  data  than  in  the  static  sensor  values.  Reasons 
for  this  included  sensor  dynamics,  cross-track  measurement  errors,  and  the  feed¬ 
back  nature  of  the  system.  The  Simulink  block  diagrams  used  to  model  the  incli¬ 
nometers  are  shown  in  Figure  19,  with  the  top  view  being  the  entire  inclinometer 
block  and  the  bottom  being  one  of  the  inclinometer  error  blocks. 


1.  Sensor  Dynamics 

All  the  sensors  used  in  the  model  (accelerometers,  rate  gyros,  and  inclinom¬ 
eters)  are  essentially  low  pass  filters;  their  ability  to  measure  desired  quantities 
diminishes  with  the  speed  with  which  these  quantities  are  changing.  With  error 
data  from  [Ref.  1]  and  bandwidth  data  from  [Ref.  2],  the  inclinometers  were  mod¬ 
eled  as  frrst-order  filters  with  a  cutoff  frequency  of  2  Hertz.  These  inclinometers 
had  a  time  constant,  or  time  for  the  output  to  reach  63%  of  the  input  value,  of  0.5 
seconds;  therefore,  any  angular  change  was  not  instantaneously  measured.  One 
way  used  to  overcome  this  difficulty  was  to  complement  the  low  frequency  incli¬ 
nometers  vidth  a  higher  frequency  device.  On  the  Bluebird  model,  this  was  accom¬ 
plished  by  integrating  the  higher  frequency  rate  gyros  and  adding  the  integrated 
output  to  the  respective  inclinometer  ouq)ut  [Ref.  3].  Both  the  unfiltered  inclinom¬ 
eter  ouqiut  and  the  filtered  inclinometer-rate  gyro  measurement  were  compared 
with  the  actual  aircraft  states  to  determine  errors  in  roll,  pitch,  and  yaw  at  each 
point  in  time.  The  resulting  errors  are  discussed  in  the  following  chapter. 
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2.  Cross  Axis  Measurement  Errmrs 


Cross  axis  error  refers  to  the  false  measurement  signal  in  one  channel  due  to 
an  actual  quantity  in  another  channel.  For  example,  a  pure  pitch  rate  will  induce 
some  small  indicsUion  in  the  roll  and  yaw  rate  gyros  as  weU.  For  the  Bluebird 
mo(tel,  the  mt^nitude  of  cross-axis  coupling  was  0.5%  for  the  accelerometers  and 
rate  gyros  [Ref.  2]. 


3.  Feedback 

In  an  open  or  closed  loop  linear  system,  statistical  analysis  of  normally  dis¬ 
tributed  errors  either  singly  or  in  combination  is  a  relatively  straightforward  {xo- 
cess.  When  one  considers  a  closed-loop  nonlinear  system  in  which  states  and 
errors  are  being  fed  back  into  the  system,  the  well-known  relationships  which  held 
for  the  linear  case  are  no  longer  valid.  To  analyze  errors  in  the  system,  one  must 
either  form  a  linearized  version  or  look  at  the  system  outputs  over  many  trials,  the 
so-called  Monte  Carlo  method.  The  first  approach,  forming  a  linearized  model  for 
analysis,  was  not  considered  because  of  the  inherent  error  introduced  in  lineariza¬ 
tion. 


VI.  RESULTS 


This  chapter  contains  results  for  the  simulation  performed  in  Chapter  5.  Monte 
Carlo  analysis  was  performed  across  the  entire  set  of  7500  data  points  for  each  run 
to  obtain  the  probability  distribution  functions  for  navigation  and  angular  errors. 
At  each  time  step  of  .01  seconds,  the  navigation  error  was  obtained  by  taking  the 
difference  between  the  inertial  position  calculated  in  the  Equations  of  Motion 
Block  and  the  position  estimate  obtained  bom  the  Navigation  Block.  Sensor  errors 
were  obtained  by  subtracting  the  outputs  of  the  sensor  models  from  the  sensor  val¬ 
ues  in  the  equations  of  motion.  The  standard  deviations  of  the  navigation  and  sen¬ 
sor  errors  are  shown  in  Table  3.  The  navigation  errors  in  Table  3  are  the  differences 
between  the  estimated  position  and  the  position  calculated  in  the  aircraft  equations 
of  motion  algorithm.  The  first  row  of  die  aircraft  and  payload  sensor  errors  is  the 
la  value  of  the  noise  fed  into  the  system.  The  second  row  is  the  standard  deviation 
calculated  from  the  set  of  data  points  obtained  using  Monte  Carlo  Simulation. 

To  validate  the  randomness  of  the  errors  introduced  into  the  model,  the  payload 
sensor  errors,  which  were  not  subject  to  any  of  the  dynamic  conditions  that  the  used 
vehicle  sensors  were,  are  also  included  in  the  statistical  analysis  of  the  collected 
data.  Assuming  a  sufficiendy  large  random  sample,  the  standard  deviation  of  the 
analyzed  payload  sensor  errors  should  match  the  standard  deviation  of  the  errors 
introduced  into  the  payload  model,  as  it  does. 
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TABLE  3:  NAVIGATION  AND  SENSOR  ERRORS 


Navigation 

(Feet) 

Aircraft 

Sensors 

(Degrees) 

Payload 

Sensors 

(EJegrees) 

X 

Y 

Z 

Roll 

Htch 

Yaw 

RoU 

Pitch 

Yaw 

Nominal 

Value 

.10 

.10 

.10 

.285 

.285 

.285 

Derived 

Value 

4.2 

6.2 

16.4 

.118 

.126 

.105 

.287 

.284 

.286 

The  nominal  and  derived  sensor  errors  were  placed  in  the  UAV  Error  Analysis 
Program  provided  by  NADC  Warminster  [Ref.  13]  and  a  comparison  was  macte  of 
the  error  in  targeting  a  fixed  site  assuming  a  one  mile  lateral  standoff  and  2000  feet 
of  vertical  separation  from  the  taiget  The  Error  Analysis  program  assumed  a  flat 
ground  plane  around  the  taiget  and  pnoject^  a  spheroid  of  one  standard  deviation 
of  errors  in  inertial  X,  Y,  and  Z  directions  from  the  UAV  onto  the  ground  plane. 
The  resulting  projection  formed  an  elliptical  curve  on  the  ground  plane.  Data  from 
the  error  analysis  program  is  presented  in  Table  4.  The  first  column  contains  the 
noise  values  used  in  the  sensor  and  payload  models.  The  next  column  contains  the 
values  obtained  from  analysis  of  the  simulation  output  data.  In  the  third  column  are 
the  sensitivities  for  each  parameter  (how  much  an  error  in  the  given  parameter  con¬ 
tributes  to  Circular  Error  Probability,  or  CEP).  The  CEP  is  the  size  of  the  radius 
placed  around  the  taiget  which  would  enclose  50%  of  the  estimates  of  target  posi¬ 
tion.  The  CEP  for  the  nonlinear  model  was  about  1 1  %  higher  than  when  consider¬ 
ing  only  the  static  sensor  errors.  In  order  to  get  the  CEP  of  the  static  model,  the 
UAV  would  have  to  get  ^roximately  10%  closer  to  the  target  at  its  present  2(XX)  ft 


51 


altitude.  From  the  sensor  accuracies  for  various  architectures  in  [Ref.  1],  tte  model 
with  the  ii^orporated  MIAG  errors  is  superior  in  accuracy  to  the  other  architec¬ 
tures. 


Table  4:  ERROR  ANALYSIS  PROGRAM 


ERRORS 

SENsmvmES 

CONTRIBUTION 

NOM 

ACT 

NOM 

ACT 

Azimuth 

.285 

.285 

32.7  m/deg 

4% 

3% 

Elevation 

.285 

.285 

111.6  m/deg 

43% 

38% 

Altitude 

45’ 

45’ 

.9  m/ft 

6% 

5% 

Roll 

1 

1.18 

79.2  m/deg 

22 

19% 

Kteh 

1 

1.26 

79.2  m/deg 

22% 

30% 

Heading 

1 

1.05 

32.7  m/deg 

4% 

4% 

North 

Error 

43.1  m 

46.1  m 

East  Error 

36.8  m 

42.2  m 

CEP 

56.7  ra 

62.5  m 

Table  4  also  shows  how  error  sensitivities  changed  with  the  analysis  of  closed 
loop  sensor  errors.  The  error  sensitivities  are  potentially  very  important  parameters 
when  consirteiing  various  UAV  architectures  for  target  acquisition  performance. 
From  the  position  estimate  errors  and  the  angular  measurement  errors  in  both  the 
aircraft  and  payload  sensors,  a  Matlab  routine,  DATERR.M  was  used  to  calculate 
position  error  at  each  point  in  time.  To  do  this,  data  for  one  run  was  Hrst  loaded 
into  memory.  At  each  time  step,  the  instantaneous  heading,  lookdown,  and  roil 
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angles  to  the  taiget  H'q.Sq,  and  d>o  were  calculated  from  taiget  position  (x^.y^z,)  and 
aircraft  true  position  (X(),yo,zo)  using 


yt-yo 


=  atan 


^t-^0 

yt-yo 


(EQ48) 


Next,  the  inertial  attitude  errors  x^^,  y^,  and  were  calculated  using  the  relationship 
relating  inertial  distances  x,  y,  and  z  to  Euler  angles  and  true  distance  to  target  R: 


X 

cosTqCOsOq 

y 

sin'E^cosO^ 

z 

sinOQCOs0Q 

(EQ49) 


For  small  angular  errors  A\|f,  A6,  and  A(t>,  we  make  the  approximation 


COS'PqCOsSq  -  cos  ( ^0  ~  cos  (®0  ~ 
sin'PQCos^Q  -  sin  (Vq  -  Ay)  cos  (Oq  -  A0) 
sinOgCosSg  -  sin  -  A^)  cos  (0q  -  A0) 


(EQ50) 


The  navigation  errors  Xq,  yg,  and  are  found  by  simply  taking  the  difference 
between  the  afrcraft  position  (xo,yo,Zo)  and  its  estimate: 
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(EQ51) 


*0 

X 

= 

yo 

— 

7 

"n 

*0 

z 

The  components  of  the  navigation  and  adimde  errors  were  then  atkkd  to  obtain  the 
entire  targeting  error.  Figure  20  is  a  plot  of  targeting  error  magnitude  for  one  run 
with  components  normalized  by  distance  to  target  The  vertical  axis  rejxesents  tar¬ 
get  error  while  the  horizontal  axis  is  time  in  huiklredths  of  seconds.  The  Probabil¬ 
ity  Density  Functions  for  navigation  error,  normalized  attitude  error,  arul 
normalized  total  error  are  shown  in  Figures  21  and  22.  These  two  graphs  are  for 
total  error  magnitude  and  do  not  consider  tnor  biases  in  any  direction,  which  were 
very  nearly  zero  for  bodi  navigati<Hi  and  angular  errors. 

An  important  consideration  when  analyzing  a  UAV  is  how  onboard  sensors  are 
integrated  to  obtain  data.  In  the  simulation  of  Bluebird,  angular  information  was 
obtained  frcnn  both  the  inclinometers  alone  and  from  relatively  slow  inclinometers 
filtered  with  larger  bandwidth  rate  gyro  information  [Ref.  3].  Table  5  presents  the 
filtered  and  unfUtered  standard  deviations  and  the  percent  improvement  for  the  fil¬ 
tered  model. 
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TUMe  5:  COMPARISON  OF  FILTERED  AND  UNnLTERED  ANGULAR  DATA 


Angle 

Measured 

UnfUtered 
Value,  la, 
degrees 

Filtered 
Value,  la, 
degrees 

Percent 

Improvement 

RoU 

1.24 

1.04 

19 

Pitch 

1.38 

1.24 

12 

Heading 

1.08 

0.36 

67 

The  laige  improvement  in  heading  accuracy  resulted  because  there  was  no 
delay  modeled  in  the  heading  sensor  [Ref.  2].  The  results  from  Table  5  are  pre¬ 
sented  in  gn4>hical  form  in  Figure  23.  For  each  axis,  the  addition  of  a  rate  sensor 
improved  angular  measurement  performance.  Figure  24  shows  the  effect  on  CEP 
of  the  filtered  and  unfiltered  sensors.  Not  only  are  the  accuracies  of  individual  sen¬ 
sors  important;  the  way  that  a  sensor  package  is  integrated  has  a  bearing  on  how 
accurately  measurements  may  be  made.  The  comparison  was  made  using 
[Ref.  12].  Standoff  distance  was  one  mile  and  vertical  separation  from  the  target 
was  2000  feet.  Other  error  inputs  to  the  program  were  used  from  the  analysis  val¬ 
ues  in  column  2  of  Table  4. 


55 


Figure  20.  Tugeting  Error  NormaUzed  By  Distance 


COMPARISON  OF  ANGULAR  ERRORS 


I  Unfilterec  □  Filtered  H  Nomina 

■  ■!■■■ 

Roll  Pitch  Yaw 

Figure  23.  Angular  Error  Comparison 


CEP  COMPARISON  FOR  DIFFERENT  SENSOR 

ERRORS 


NOMINAL  UNFILTERED  FILTERED 


Figure  24.  CEP  Comparison 


vn.  CONCLUSIONS  AND  RECOMMENDATIONS 


Wl^n  comparing  different  dynamic  systems  with  regard  to  how  well  they  per¬ 
form  a  certain  mission  given  onboard  sensors  and  their  accuracies,  it  is  important  to 
look  beyond  how  sensors  perform  independently  in  a  static  situation  and  to  con¬ 
sider  how  they  perform  in  a  dynamic  environment  together  with  the  rest  of  a  vehi¬ 
cle’s  sensor  package.  The  key  to  this  analysis  is  the  formulation  of  a  high  fidelity 
model  containing  all  sensor  and  vehicle  motion  dynamics.  Because  such  a  model 
will  necessarily  incorporate  many  nonlineaiities,  the  most  general  approach  in 
accurately  analyzing  system  performance  is  Monte  Carlo  simulation. 

In  the  UAV  model  and  targeting  simulation  considered  in  this  thesis,  it  was 
shown  that  static  sensor  accuracies  differed  greatly  from  measured  data  obtained 
from  the  sensors  during  simulation.  Also,  it  was  important  to  look  at  what  sensors 
were  being  used  to  provide  data.  For  example,  whether  low  frequency  sensors 
were  complemented  with  higher  frequency  ones  to  provide  the  vehicle  with  better 
state  measurements  during  maneuvering  flight. 

The  Naval  Postgraduate  School  Avionics  Laboratory  has  adequate  facilities 
now  to  create  high  fidelity  models  and  run  them  in  realistic  simulations.  As  the 
models  created  by  successive  graduate  students  with  little  {n-ogramming  experience 
become  more  and  more  complex,  time  constraints  require  that  models  created  in  a 
user-friendly  graphical  environment  be  able  to  be  run  in  as  optimal  an  environment 
as  possible.  The  Aeronautics  Department  recently  acquired  MATRIXx.  a  package 
which  allows  construction  of  models  in  a  graphical  environment  similar  to  that  of 
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Simuiink  but  with  an  autocode  feature  which  creates  executable  C  code  from  the 
graphical  model  to  allow  for  faster  simulation. 

In  the  near  term,  the  model  used  in  this  paper  presents  a  good  teaching  aid  for 
both  avionics  system  design  and  integration  and  LQR  control  design.  For  follow- 
on  thesis  students,  the  modularity  of  the  UAV  model  allows  for  the  modification  or 
complete  reconstruction  of  any  or  all  of  the  subsystems.  This  provides  an  ideal 
environment  for  future  studies  in  advanced  controller  design  and  other  avionics 
topics.  Suggested  topics  include: 

1 .  The  implementation  of  a  fault  detection  system  on  the  model. 

2.  The  design  of  a  controller  based  solely  on  inertial  states. 

3.  The  redesign  of  the  sensor  package. 

The  increasing  cost  of  manned  aircraft  and  the  continuing  miniaturization  of 
avionics  components  will  make  Unmanned  Aerial  Vehicles  more  and  more  attrac¬ 
tive  as  an  alternative  to  airplanes  and  helicopters  for  many  missions.  It  is  incum¬ 
bent  that  the  government  maintain  adequate  facilities  to  realistically  evaluate  future 
UAV  designs.  The  Naval  Postgraduate  School  Avionics  Laboratory,  with  its  mix  of 
industry-experienced  professors  and  fleet  aviators,  provides  an  excellent  environ¬ 
ment  for  such  work. 
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APPENDIX  A:  MATLAB  ROUTINES 


A.  GUIDANCE  ROUTINE 
%  Input  vector  is  psi,  time 
%  ou4)uts  incremental  px,  py  for  this  iteration 
%  also  ouQmts  place  holder  for  guidance  mechanism 
%define  function 

function  out  =  psicom2(u,pxc,pyc,told,ss) 

%define  waypoints  in  inertial  X,Y 
pxcv=[1000  6000  100000] 
pycv=[1000  3000  30000] 

%  inputs  to  block  are  heading,  time,  X  and  Y  position,  velocity,  and  current 

%  waypoint 

u=u' 

psi=su(l) 

px=u(3) 

py=u(4) 

speeds:u(5) 

mew=u(2) 

%  start  using  waypoint  at  0.2  seconds  of  flight  time 
if  tnew>.2 
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ss3u(6) 

elsesssl 

end 

%  define  commanded  waypoint 

pxc=spxcv(ss) 

pyc=pyc(ss) 

yy=ss; 

%  step  to  next  waypoint  if  within  2  seconds  of  current 
peiTOr=ssqrt((pxc-px)^2+(pyc-py)^2) 
if  perror  <=  150 
yysyy+1 
end 

%  calculate  ground  track  heading  to  waypoint 

psic=-atan2(pyc-py,pxc-px) 

out=:[psic  yy]' 


B.  LQR  DESIGN  ROUTINE 

%  load  a,b  matrices 
[Ai,Bi,CiJ)i]=:linmod('newsynt2'); 

%  B2  is  the  input  matrix  across  the  actuators 
B2=Bi(:,5:8); 

%B1  is  the  command  input  matrix 
Bl=:Bi(:,l:4) 
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%  Cl  is  the  ou^t  matrix  for  the  three  commanded  states-airspeed,  attitude, 
%  and  ground  track  heading 

C1=[0000  1  zeros(l,12);zeros(l,ll)  1  00000;zeros(l,12)  1  0000]; 
pause 

%  solve  algebraic  Ricatti  equation  to  get  state  feedback  gains 
p=aie(Ai.B2*inv(R)*B2',Ci'*Q*Ci); 

%  kstate  is  state  feedback  gains 

%  kerror  is  command  error  gain  matrix 

k=inv(R)*B2’*p; 

kstate=k(:,l:13) 

kenor=k(:,14:17) 

%  get  closed  loop  eigenvalues 
eig(Ai-B2*k); 

%  get  damping  of  closed  loop  eigenvalues 

damp(ans) 

pause 

%  set  bode  frequency  range  from  .01  to  100  rad/s 

w=logspace(-2,2); 

wl=w'; 

%  solve  frequency  response  for  velocity/velocity  commanded 
[m  1  ,pl  ]=bode(  Ai-B2*k,B  1  ,C  1  ,zeros(3,4),2,  w); 

%  convert  to  dB 
ml=20*logl0(ml); 
mmls:[wl  ml]; 
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save  datal  mml  /ascii 


%  size  plot  for  convenient  size  when  creating  postscript  file  for  thesis 
subplot(1.3,1.3,1.2) 

%  full-scale  bode  magnitude  plot  for  velocity/velocity  commanded 

semilogx(w,ml) 

grid 

xlabeK'FREQUENCY  (RAD/S)*) 
ylabelCdB  MAGNITUDE) 
titleCVelocity  to  Velocity  G>mmanded') 
pause 

%  create  bode  magnitude  plot  for  attitude/attitude  commanded 

[m2,p2]=bode(  Ai-B2*k,B  1  ,C  1  ,zeros(3,4),3,w); 

m2=20*logl0(m2); 

mm2s[wl  m2]; 

save  data2  mm2  /ascii 

subplot(1.3,1.3,1.2) 

semilogx(w,m2) 

grid 

xlabeK’FREQUENCY  (RAD/S)') 
ylabelCdB  MAGNITUDE  ) 
tideC  Attitude  To  Attitude  Commanded') 
pause 
% 

%  create  bode  plot  for  ground  track/ground  track  commanded 
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{m3,p3]sbode(  Ai'B2*k,B  1  ,C1  ,zeros(3,4),4,  w); 

m3*20*log  10(m3); 

mm3=[wl  m3]; 

save  data3  mmS  /ascii 

subplot(1.3,1.3,1.2) 

semilogx(w,m3) 

grid 

xlabeK'FREQUENCY  (RAD/S)) 

ylabelCdB  MAGNITUDE) 

titieCGround  Track  to  Ground  Track  Commanded') 

pause 

%  create  bode  plot  for  elevator/elevator  commanded 

[m4,p4]asbodc(  Ai-B2*k,B2,k,zeros(4,4),  I  ♦  w); 

m4=20*log  10(m4); 

mm4s[wl  m4]; 

save  data4  mm4  /ascii 

subplot(1.3,1.3,1.2) 

semilogx(w,m4) 

grid 

xlabelCFREQUENCY  (RAD/S)') 
ylabelCdB  MAGNITUDE  ) 
tideCElevator  to  Elevator  Commanded') 
pause 

%  create  bode  plot  for  rudderAudder  commanded 
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[mS,pS]sbode(Ai-B2’''k,B2,k,zeros(4,4),2,w); 

m5=20*log  10(m5); 

inm5=[wl  m5]; 

save  dataS  mm5  /ascii 

subplot(1.3,1.3,1.2) 

semilogx(w,mS) 

grid 

xlabeK’FREQUENCY  (RAD/S)’) 
ylabeK’dB  MAGNITUDE') 
tideCRudder  to  Rudder  Commanded') 
pause 

%  create  bode  plot  for  aileron/aileron  commanded 

[m6,p61s=bode(Ai-B2*k,B2,k,zeros(4,4),3»w); 

m6=20*logl0(m6); 

mm6=[wl  m6]; 

save  data6  mm6  /ascii 

subplot(1.3,1.3,1.2) 

semilogx(w4n6) 

grid 

xlabel('Frequency  (Rad/S)') 
ylabeK’dB  MAGNITUDE’) 
tideC  Aileron  to  Aileron  Commanded') 
pause 

%  bode  plot  for  throtde/throtde  commanded 
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[m7,p7J=bodc(Ai-B2*k,B2,k.zeros(4,4X4,w); 

m75s20*Iog  10(m7); 

iiun7=[wl  m7]; 

save  data7  min7  /ascii 

subplot(1.3,1.3,1.2) 

semilogx(w4n7) 

grid 

xiabelCFREQUENCY  (RAD/S)') 
ylabelCdB  MAGNITUDE) 
tide('Thiottle  to  Throttle  Commanded’) 
pause 

%  create  time  vector  for  time  response  plots 

t=[0:.5:501; 

tlsf; 

%  step  response  from  velocity  command  to  velocity 

m8=step(Ai-B2*k,B  1  ,C1  ,zeros(3,4),2,t); 

mm8=[tl  m8]; 

save  dataS  mm8  /ascii 

subplot(1.3,1.3,1.2) 

plot(t,m8) 

grid 

titleCResponse  to  a  Step  Velocity  Command’) 
xlabeK’TIME  (SECONDS)’) 
ylabeK’RESPONSE') 
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pause 

%  step  response  from  attitude  command  to  attitude 

m9=step(Ai-B2*k,B  1  ,C1  ,zeros(3,4),3.t); 

nun9s:[tl  m9]; 

save  data9  mm9  /ascii 

sub{dot(1.3,1.3,1.2) 

plot(t,m9) 

grid 

title('Response  to  a  Step  Attitude  Command') 
xlabelCTIME  (SECONDS)’) 
ylabel(’RESPONSE) 
pause 

%  ground  track  heading  step  response  to  ground  track  input 

ml0s=step(Ai-B2*k,B  1  ,C  1  ,zeros(3,4),4,t); 

mmlOs[tl  mlO]; 

save  datalO  mmlO  /ascii 

subplot(1.3,1.3,1.2) 

plot(t,mlO) 

grid 

titleCResponse  to  a  Step  Heading  Command') 
xlabeKTlME  (SECONDS)') 
ylabeK'RESPONSE) 
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C.  TARGETING  ERROR  CALCULATION 
%  error.m 

%  This  Matlab  routine  takes  saved  bliKbird  navigational  data  and  computes 
%a  targeting 

%  error  at  .01  second  intervals 
% 

%  Load  data  from  run 
loaddl 

%  Target  is  at  (6000,0,0)  feet  in  (X, Y,Z)  universal  coordinate  system 

xt=6000 

ytsO 

zt=0 

%  Make  error  computation  at  each  step 
fori=l;7500 
%  Aircraft  true  position 
x0=P(i,10); 
y0=P(i,ll); 

20=P(i,12); 

%  True  heading  to  target 
psi0=-atan2(-y0,xt-x0); 

%  True  lookdown  angle  to  target 
%  Assuming  2000  ft  altitude 
thetaO=:>atan2(zO,xt-xO); 

%  True  roll  angle  to  target 
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th- 

\ 

I^Os:ataii2(zO,-yO); 

%  True  range  to  target 
R(i)=sqiK(3(X)0-x0)^2+y0^2+z0^2); 

% 

%  Get  instantaneous  angular  errors  from  file 
dirfiisstates(i,7); 
dtheta=states(i,8); 
dpsi=states(i,9); 

% 

%  Calculate  inertial  X,Y,Z  errors  from  attitude  sensors 
xerratt(i)=R(i)'''(cos(psiO)*cos(thetaO)'Cos(psiO+dpsi)'^cos(theta(>fdtheta)); 
yeiTatt(i)=:R(i)*(sin(psiO)*cos(phiO)-sin(psiOfdpsi)*cos(phiOK!phi)); 
zeiTatt(i)=R(i)*(sin(phiO)*cos(thetaO>sin(phiO+dphi)*cos(thetaO+dtheta)); 
% 

%  Calculate  inertial  X,Y,Z  errors  from  navigation 
xemiav(i)=P(i,  1 0)-Phat(i,  1 ); 
yermav(i)=P(i,  1 1  )-Phat(i,2); 
zermav(i)=P(i,  1 2)-Phat(i,3); 

% 

%  Sum  errors  in  X,YZ 
xerr(i)=xermav(i)+xerratt(i); 
yerr(i)5syermav(i)+yerratt(i); 
zerr(i)=zermav(i)+zerratt(i); 

%  Total  error 
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iientot(i>=sqrt(xeiTnav(i)^2+yermav(i)^2+2ermav(i)^2); 

errtot(i>5sqit(xen(iy'2+yerr(iy'2+zerr(i)^2); 

miieiT(i)=:errtot(i)/R(i); 

end 


APPENDIX  B.  MATHEMATICAL  PROPERTIES 


A.  LINEARIZED  BLUEBIRD  PLANT  MATRICES 
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B.  GROUND  TRACK  HEADING  DERIVATION 

Consider  an  aircraft  flying  in  straight  and  level  flight  with  heading  y.  Now  add 
a  sideslip  angle  Because  ^  is  in  body  coordinates,  it  must  be  converted  to  inertial 
coordinates  to  add  its  effect  to  heading.  For  a  cruise  condition,  this  is  accomplished 
by  rotating  ^  into  the  inertial  X-Y  plane  by  the  climb  angle  y. 
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Vcosy  =  -p. 


(EQ47) 


A  positive  sideslip  angle  p  corresponds  to  a  left  crab  angle,  while  heading  angle  is 
positive  to  the  right. 

The  roll  angle  (|>  also  contributes  to  ground  track  heading.  Since  (j)  is  measured 
in  body  coordinates  as  well,  it  must  also  be  rotated  to  inertial  coordinates  by  the 
climb  angle.  It  contributes  to  \|r  with  angle  of  attack.  A  positive  roll  angle  corre¬ 
sponds  to  a  right  wing  down  condition. 


\|rcosY  =  (]>sina  (eq48) 

Heading  angle  increases  to  the  right  when  viewing  the  inertial  X-Y  plane  from 
above.  Since  sideslip  and  roll  angle  are  used  to  counteract  drift,  they  contribute  to 
ground  track  heading  positively  and  negatively  respectively.  Therefore,  the  rela¬ 
tionship  used  to  obtain  ground  track  heading  becomes: 


¥+ 


p-<j)sina 

cosy 


(EQ  49) 


C.  SIMUJ  ATION  STARTUP  PROCEDURES 

1.  Logon 

2.  In  UNIX  Cmnmand  Window,  type  "xhost  +" 

3.  Set  environment  for  the  machine  in  use: 

"setenv  DISPLAY  13L120.149.xx"  where  xx  is  a  machine-specific  code. 

4.  type  "matlab" 
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5.  Change  to  the  £^ropriate  directory. 

6.  type  "birdset".  This  is  a  .M  file  which  loads  startup  variables  into  the  work¬ 
space  and  calls  the  Simulink  program  BIRO.NLS.M. 

7.  From  the  BIRD_NL5.M  pulldown  menu,  select  "Simulation",  then  "Start". 

8.  Variables  saved  to  the  workspace  include  P  (actual  positions  and  states),  Phat 
(estimated  positions),  and  states(differences  between  estimated  and  actual  states). 

9.  Simulation  runs  on  Hornet  (a  Sparc  10  workstation)  at  a  rate  of  one  second  of 
real  time  for  every  five  minutes  of  real  time. 
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